function distortedPoints = UndistortedPoints(points, intrinsicMatrix, radialDistortion, tangentialDistortion)   
   %#codegen
    %视觉十四讲 pdf135页
    % unpack the intrinisc matrix   
    cx=intrinsicMatrix(1,3);
    cy = intrinsicMatrix(2,3);
    fx =intrinsicMatrix(1,1);
    fy = intrinsicMatrix(2, 2);
    skew = intrinsicMatrix(1, 2);

    % center the points
    center = [cx, cy];
    centeredPoints = points-center;
    
    % normalize the points
    yNorm = centeredPoints(:, 2, :) ./ fy;
    xNorm = centeredPoints(:, 1, :) ./ fx;
    
    % compute radial distortion
    r2 = xNorm .^ 2 + yNorm .^ 2;
    r4 = r2 .* r2;
    r6 = r2 .* r4;
    
    k = zeros(1, 3, 'like', radialDistortion);
    k(1:2) = radialDistortion(1:2);
    if numel(radialDistortion) < 3
        k(3) = 0;
    else
        k(3) = radialDistortion(3);
    end
    
    alpha = k(1) * r2 + k(2) * r4 + k(3) * r6;
    
    % compute tangential distortion
    p = tangentialDistortion;
    xyProduct = xNorm .* yNorm;
    dxTangential = 2 * p(1) * xyProduct + p(2) * (r2 + 2 * xNorm .^ 2);
    dyTangential = p(1) * (r2 + 2 * yNorm .^ 2) + 2 * p(2) * xyProduct;
    
    % apply the distortion to the points
    Xcnp=[xNorm, yNorm];                         % Xcn'
    normalizedPoints=cell(1,20);                   % Xcn
    delta_distortedNormalizedPoints=cell(1,20);    % delta_Xcn
    for m=0:20
        if m==0
            normalizedPoints{1,1}=Xcnp;
            delta_distortedNormalizedPoints{1,1}(1,1:2)=0;
        end

        if m>0
            normalizedPoints{1,m+1}=Xcnp-delta_distortedNormalizedPoints{1,m};
            xNorm= normalizedPoints{1,m+1}(:,1);
            yNorm= normalizedPoints{1,m+1}(:,2);
            xyProduct = xNorm .* yNorm;
            r2 = xNorm .^ 2 + yNorm .^ 2;
            r4 = r2 .* r2;
            r6 = r2 .* r4;
            alpha = k(1) * r2 + k(2) * r4 + k(3) * r6;

           dxTangential = 2 * p(1) * xyProduct + p(2) * (r2 + 2 * xNorm .^ 2);
           dyTangential = p(1) * (r2 + 2 * yNorm .^ 2) + 2 * p(2) * xyProduct;
           delta_distortedNormalizedPoints{1,m+1} =  normalizedPoints{1,m+1}.* [alpha, alpha] + ...
            [dxTangential, dyTangential];  
        end   
    end

    % convert back to pixels
    distortedPointsX = normalizedPoints{21}(:, 1, :) * fx  + skew * normalizedPoints{21}(:, 2, :);
    distortedPointsY = normalizedPoints{21}(:, 2, :) * fy;  
    distortedPoints = [distortedPointsX, distortedPointsY]+center;

end